Hadron spectra and elliptic flow for 200 A GeV Au+Au collisions from viscous 
hydrodynamics coupled to a Boltzmann cascade 
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It is shown that the recently developed hybrid code VISHNU, which couples a relativistic viscous 
fluid dynamical description of the quark-gluon plasma (QGP) with a microscopic Boltzmann cascade 
for the late hadronic rescattering stage, yields an excellent description of charged and identified 
hadron spectra and elliptic flow measured in 200 A GeV Au+Au collisions at the Relativistic Heavy- 
Ion Collider (RHIC). Using initial conditions that incorporate event-by-event fluctuations in the 
initial shape and orientation of the collision fireball and values n/s for the specific shear viscosity of 
the quark-gluon plasma that were recently extracted from the measured centrality dependence of the 
eccentricity-scaled, pr-integrated charged hadron elliptic flow W2, c h/e, we obtain universally good 
agreement between theory and experiment for the pr-spectra and differential elliptic flow V2(pr) for 
both pions and protons at all collision centralities. 

PACS numbers: 25.75.-q, 12.38. Mh, 25.75.Ld, 24.10.Nz 
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I. INTRODUCTION 

In a recent article [1] we extracted the shear viscos- 
ity to entropy density ratio (?7/s)qgp of the quark-gluon 
plasma (QGP) created in heavy-ion collisions at RHIC 
by comparing experimental data for the eccentricity- 
scaled elliptic flow t>2 1 £ with calculations performed with 
VISHNU [2], a hybrid model that describes the QGP stage 
of the expansion of the collision fireball macroscopically 
with viscous hydrodynamics (in which {t]/s)qqp enters 
as an input parameter) but switches to a microscopic de- 
scription in the late hadronic phase where we solve the 
Boltzmann equation with UrQMD [3]. An important step 
in this analysis was to ensure that, as we compared the- 
oretical curves for different (j]/s)qgp values with the ex- 
perimental data in order to find the value preferred by 
Nature, we maintained a good description of the total 
charged hadron multiplicity and the hadron transverse 
momentum spectra as a function of collision centrality. 
That this was indeed achieved was announced in Rcf. [1] 
and will be documented in this companion article. Wc 
then proceed to demonstrate that, with the QGP shear 
viscosity extracted in [1], VISHNU provides a good descrip- 
tion of all single-particle aspects of soft hadron produc- 
tion in 200 A GeV Au+Au collisions for which accurate 
measurements exist, over the entire range of collision cen- 
tralities. 



II. METHODOLOGY 

The various components of the viscous hydrodynam- 
ics+Boltzmann hybrid code VISHNU have been described 
in Rcf. [2] (see also [3-5]) to which we refer the reader 
interested in technical details. For the QGP fluid we ap- 
proximate rj/s in the temperature range T c < T < 2T C by 



a constant [6]. We switch from a hydrodynamic descrip- 
tion of the QGP to the microscopic hadronic rescattering 
code UrQMD at temperature T sw = 165 MeV, adjusted to 
reproduce the chemical freeze-out temperature measured 
in RHIC collisions [7] ; as shown in [2] , this is at the same 
time the highest T for which we have a valid microscopic 
description and the lowest T for which the macroscopic 
hydrodynamic approach can be trusted. 

As shown in Ref. [1, 8], the QGP shear viscosity 
( 7 ?/ s )qgp extracted from the experimentally measured 
elliptic flow depends on the initial fireball eccentricity 
£part = where x and y label the coordinates along 

the short and long major axes of the fireball in the plane 
transverse to the beam direction. (This definition of x, 
together with the beam direction z, define the "partici- 
pant plane", reflected in the subscript.) With presently 
available tools this initial eccentricity cannot be directly 
measured, and theoretically we have limited control over 
it. We here use initial entropy density profiles from two 
popular geometric models for the initial particle produc- 
tion in high-energy heavy-ion collisions, the Monte Carlo 
Glauber model (MC-Glaubcr [9]), in a version [10] that 
uses finite size nucleons, and the Monte Carlo fKLN (MC- 
KLN) model [10-12]. These models give initial eccen- 
tricities that differ (depending on centrality) by up to 
25% which we hope to cover the physically reasonable 
range of uncertainty. We showed in [1] that this uncer- 
tainty in the initial eccentricity completely dominates the 
present error range in the phenomenological extraction of 
( 7 ?/ s )qgp 1 an d that future improvements in the accuracy 
of the experimentally extracted value of (?7/s)qgp can- 
not be achieved without obtaining better (experimental 
and/or theoretical) control over the initial fireball eccen- 
tricity. 

Due to the finite number of nucleons colliding with 
each other in a heavy-ion collision, the initial eccentric- 
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FIG. 1: (Color online) pr-spectra of pions (left) and protons (right) for 200 A GeV Au+Au collisions of different centralities as 
indicated. Data from the STAR (x, [15-17]) and PHENIX (+, [18]) experiments are compared with VISHNU calculations using 
MC-Glauber (dashed lines) or MC-KLN initial conditions (solid lines) and different values r\/s for the QGP shear viscosity as 
indicated. Different rj/s values are associated with different starting times to for the hydrodynamic evolution as discussed in 
the text. The STAR and PHENIX proton data shown in the right column are feeddown-corrected by removing protons from 
weak hyperon decays [15, 18]. Where necessary, PHENIX yields from neighboring narrower centrality bins were averaged to 
obtain data in the wider centrality bins used by the STAR Collaboration. 



ity of the density of secondary particles produced in these 
collisions fluctuates from event to event, as does the ori- 
entation of its major and minor axes relative to the re- 
action plane [13] (defined by the directions of the impact 
parameter and the beam). To account for these event-by- 
event fluctuations on average, we use a Monte Carlo sam- 
pling procedure to generate from the Glauber and fKLN 
models a large number of initial entropy density distribu- 
tions whose shape and orientation fluctuate from event to 
event, recenter and rotate each distribution around the 
beam direction such that its short major axis x aligns 



with the direction of the impact parameter b, sort them 
into centrality bins by N pax t (the number of wounded nu- 
clcons), and then superimpose the distributions to obtain 
a smooth average density that has the correct average 
eccentricity for collisions in this centrality class. 1 The 



Strictly speaking, this procedure yields £ par t = / tt~^t wnere 
(. . .)s denotes the expectation value taken with the averaged en- 
tropy density obtained by superimposing many recentered and 
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TABLE I: Mean (e pa rt ~ (epart)) and rms ^e{2} = y (£p art ) j participant eccentricities, as well as (£p art ) 1//a , with a computed 

from the event-plane resolution R as described in [21, 28], for the STAR [19] and PHENIX [25] Au+Au experiments at different 
collision centralities. For centralities between 10 and 50%, PHENIX [25] published data in 5% centrality increments. For the 
purpose of comparing the PHENIX and STAR data we combined the PHENIX data [25] from neighboring centrality bins, by 
averaging the event-plane resolutions R and the exponents a corresponding to the two sub-bins when computing {£p art ) 1//a for 
the larger combined bin. All eccentricities are calculated with the entropy density as weight. 



centrality 


model 


£part 


\J ( £ part) 


R [19] 


STAR 
a [21, 28] 


/ a \l/a 
\ c part / 


R [25] 


PHENIX 
a [21, 28] 


i a U/a 
\ ° p art / 


0-5% 


MC-Glauber 

IV J. \y _L VJ_j1 > 


0.089 
0.097 


0.101 
0.109 


0.61 


1.52 


0.095 
0.103 


0.51 


1.66 


0.097 
0.105 


5-10% 


MC-Glauber 
MC-KLN 


0.139 
0.172 


0.153 
0.183 


0.735 


1.31 


0.144 
0.175 


0.63 


1.49 


0.146 
0.178 


10-20% 


MC-Glauber 
MC-KLN 


0.215 
0.265 


0.230 
0.277 


0.816 


1.18 


0.218 
0.267 


0.720 


1.33 


0.220 
0.269 


20-30% 


MC-Glauber 
MC-KLN 


0.299 
0.360 


0.311 
0.372 


0.843 


1.14 


0.298 
0.362 


0.743 


1.30 


0.301 
0.364 


30-40% 


MC-Glauber 
MC-KLN 


0.361 
0.434 


0.378 
0.447 


0.825 


1.16 


0.364 
0.436 


0.704 


1.36 


0.367 
0.439 


40-50% 


MC-Glauber 
MC-KLN 


0.414 
0.493 


0.433 
0.509 


0.771 


1.25 


0.419 
0.497 


0.617 


1.50 


0.424 
0.501 


50-60% 


MC-Glauber 
MC-KLN 


0.458 
0.541 


0.481 
0.561 


0.677 


1.41 


0.468 
0.549 


0.489 


1.69 


0.475 
0.555 


60-70% 


MC-Glauber 
MC-KLN 


0.497 
0.581 


0.523 
0.606 


0.549 


1.61 


0.513 
0.597 








70-80% 


MC-Glauber 
MC-KLN 


0.528 
0.621 


0.560 
0.650 


0.412 


1.78 


0.554 
0.645 









elliptic flow resulting from the VISHNU evolution of this 
initial profile is interpreted as the event-average (x>2) for 
the selected centrality class. 

The ensemble-averaged initial entropy density is nor- 
malized such that, after evolution with VISHNU, it repro- 
duces the measured final charged hadron rapidity density 
dNch_/dy in the most central collisions [17]; due to viscous 
entropy production this is an iterative process, requir- 
ing two or three iterations. After normalization in cen- 
tral collisions, the centrality dependence of the initial en- 
tropy production is taken directly from the model (MC- 
Glauber or MC-KLN); for the MC-Glauber model we fol- 
low [10, 14] and assume a two-component (soft+hard) 
model with a small hard fraction (5 = 0.14 [10]) for the 
entropy production. In [10] this fraction was fixed within 
a hydro+Boltzmann hybrid approach using ideal fluid 
dynamics for the QGP; taking the same fraction in our 



rotated Monte Carlo events; this is not identical with, but numer- 
ically very close to the ensemble- averaged participant eccentricity 
(epart) where for each event £ par t is computed as the analogous 
expectation value taken with the entropy density of that event. 



viscous hydro+Boltzmann code ignores the centrality de- 
pendence of viscous heating. We have checked that its 
effects on the centrality dependence of the final dN^/dy 
are negligible relative to experimental uncertainties. 

III. RESULTS: COMPARISON OF SPECTRA 
AND v 2 TO DATA 

Let us now begin discussing our results. All experimen- 
tal data and theoretical calculations arc for 200 A GeV 
Au+Au collisions. Figure 1 shows pion and proton trans- 
verse momentum spectra from the STAR [15-17] and 
PHENIX [18] Collaborations for the whole range of col- 
lision centralities, separated by multiplicative factors of 
10 for clarity. The lines show VISHNU calculations for 
different values of (tj/s)qqp in the QGP phase, using 
cither MC-KLN (solid) or MC-Glauber (dashed) initial 
conditions. When changing (??/s)qgp we have to (i) 
renormalize the initial entropy density profiles by a (b- 
independent) constant factor to account for the change 
in viscous entropy production, and (ii) adjust tq to ac- 
count for the additional radial acceleration caused by the 
transverse shear pressure gradients. The latter increase 
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FIG. 2: (Color online) Eccentricity-scaled elliptic flow vije as function of pr for charged hadrons from 200 A GeV Au+Au 
collisions at different centralities. The experimental data (solid symbols) are U2{EP} measurements from the STAR experiment 
[19], scaled by {£p art ) 1 ' /c * from the Monte Carlo Glauber model (left column) and the MC-KLN model (right column), respec- 
tively, q depends on the event-plane resolution R and varies from one centrality bin to the next (see Table I). Dashed and 
solid lines with open symbols are results from VISHNU for two different values of (?7/s)qgp (0.08 and 0.16 for the MC-Glauber 
calculations, 0.16 and 0.24 for the MC-KLN calculations). The theoretical lines show the ratio (i>2)/£ P art where (...) denotes 
an average over events, and e pa rt is the eccentricity of the smooth average initial entropy density. Different symbols denote 
different collision centralities as indicated. 



with (?7/s)qgp, leading to more radial flow and flatter 
PT-spectra unless we compensate by increasing the start- 
ing time for the hydrodynamic evolution accordingly. 
The curves shown in Fig. 1 correspond to the follow- 
ing parameter pairs (r)/s, ct ): (0,0. 4fm), (0.08, 0.6 fm), 
(0.16, 0.9 fm), and (0.24, 1.2 fm). We stop at ?7/s = 0.24 
since we will see that larger QGP shear viscosities are 
excluded by the elliptic flow data. 

Except for very peripheral collisions, the different lines 
in Fig. 1 overlap almost perfectly and thus arc hard to dis- 
tinguish optically. This is intentional since it shows the 
approximate equivalence of the different parameter pairs 
as far as the quality of the theoretical description of the 
measured pr-spectra goes. Differences between theory 



and data are generally less than between data sets from 
the different experiments. We note that the theoretical 
proton spectra are uniformly about 50% larger than the 
PHEN1X data but agree nicely with their slope; their nor- 
malization agrees somewhat better with the STAR data. 
Due to limited event statistics, VISHNU does not include 
protons from weak decays; in Fig. lb we therefore com- 
pare with experimental data that have been corrected 
to eliminate feeddown protons. However, the feeddown 
correction methods used by PHENIX and STAR differ 
[15, 18], and systematic uncertainties arising from the 
feeddown correction are large. Keeping the differences 
between the experimental data sets in mind, VISHNU pro- 
vides a very acceptable compromise description. We do 
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FIG. 3: (Color online) Same as Fig. 2, but with u 2 {EP} data from the PHENIX Collaboration [25] instead of STAR data. See 
Table I for the effective eccentricities (£p art } 1//Q! used for each centrality bin. 



note in passing that in the most peripheral bins a vis- 
cous treatment of the QGP appears to work better than 
treating it as an ideal fluid; assuming zero viscosity for 
the QGP gives too little radial flow and results in 
spectra for both pions and protons that are slightly too 
steep. 

Figure 2 shows the differential elliptic flow V2{pt) for 
charged hadrons from Au+Au collisions at different cen- 
tralizes. Experimental data from the STAR Collabora- 
tion obtained with the event-plane method [19] are com- 
pared with VISHNU calculations for different QGP shear 
viscosities. ^{EP} receives positive contributions from 
event-by-event flow fluctuations and non-flow effects [20] . 
The latter can be minimized by trying to decorrelate the 
determination of the event plane from the measurement 
of i)2, e.g. by employing a large rapidity gap between 
these measurements. Fluctuation effects can not be elim- 
inated from the measurement, but both non-flow and 
fluctuations can be corrected for [21]. In [1] we used 



such corrected data to extract (t7/s)qgp; in Figure 2 we 
show the uncorrected ^{EP} data directly as measured. 
To account for the fluctuation contribution we normalize 
them by (£p ar t) 1/,a [21, 22] where the exponent a depends 
on the experimental event-plane resolution R and on de- 
tails of the V2 extraction method [21]. In Table I we have 
summarized for each centrality bin shown in Figs. 2-4 
the event-plane resolution factors R for the STAR and 
PHENIX experiments, the corresponding a values ob- 
tained from the procedure described in Ref. [21], as well 
as the corresponding values for (e pa rt), (Epart) 1 ^ 2 ' an d 
(£part) 1/Q - In Fig- 2 we compare the experimental ratio 
y2{EP}/(£p art ) 1 / Q with the theoretically calculated ratio 
(i>2)/£ P art- This is the correct comparison if v<i ~ £ pa rt 
event by event, as suggested by hydrodynamic simula- 
tions [23] (see, however, [24]). 

Figure 2 demonstrates excellent agreement between 
VISHNU and the experimental data, over the entire 
range of centralities except for the two most periph- 



6 



eral bins, if we use (?7/s)qgp = 0.08 for MC-Glauber 
(left column) and (t7/s)q GP = 0.16 for MC-KLN initial 
conditions (right column). 2 These values are a little 
smaller than, but consistent with the corresponding val- 
ues (t?/s)qgp — 0.1 for MC-Glauber and (t]/s)qgp — 0.2 
for MC-KLN that were extracted in [1] from the px- 
integrated, non-flow and fluctuation corrected charged 
hadron V2- Small non-flow effects in the ^{EP} data 
shown here, shifting them slightly upward, may account 
for this difference. 

To check this possibility, we show in Fig. 3 the same 
comparison with PHENIX data for -^{EP} [25] where the 
event plane was determined with counters several units 
of rapidity away from the central region where vi was 
measured. The PHENIX data should therefore be less 
affected by non-flow effects than the STAR data. For 
centralitics > 10% the agreement between VISHNU calcu- 
lations and the PHENIX data is equally good as in Fig. 2 
for the STAR data, with the same values for (t]/s)qgp. 
In the two most central bins, 5—10% and 0—5% re- 
spectively, the agreement deteriorates significantly, with 
the PHENIX data pointing counterintuitively to larger 
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FIG. 4: Comparison of the STAR [19] and PHENIX v 2 {EP} 
data [25] used in Figs. 2 and 3. Both data sets are scaled by 
the effective eccentricity e = {£part) 1//Q , where a depends on 
the event-plane resolution R [21] and thus varies with central- 
ity and from experiment to experiment (see Table I). 



( 7 ?/ s )qgp values for central collisions than in the other 
centrality bins. Figure 4 shows that this results from a 
disagreement between the two data sets in near-central 
collisions: while the two data sets overlap excellently for 
centralitics > 20%, they increasingly diverge at small 
centralities, with a 30% difference between STAR an 
PHENIX in the — 5% centrality bin. It has been pointed 
out that the excess of the STAR over the PHENIX data 
is uniform in px and could be explained by a 2% shift in 
the centrality definitions between the experiments [26]. 
Where such a shift could arise from and which of the two 
definitions needs to be corrected is presently under study 
[27]. Wc conclude from Figs. 2-4 that (i) non-flow ef- 
fects seem to be similar and likely small in both STAR 
and PHENIX v 2 {EP} data for centralities between 20% 
and 60%, (ii) if the STAR centrality definition is correct 
we have excellent agreement between VISHNU and the 
experimental charged hadron elliptic flow V2(pt) at all 
centralities, with (?7/s)qgp = 0.08 for MC-Glauber and 
(v/ s )qgp = 0.16 for MC-KLN initial conditions, and (iii) 
if the PHENIX centrality definition is correct, this uni- 
form agreement is broken in the most central collisions 
for which the PHENIX data appear to require larger ef- 
fective (?7/s)qgp values than at larger centralitics. 

Figure 5 demonstrates that the agreement of VISHNU 
using the phenomcnologically extracted (??/s)qgp values 
from Ref. [1] with the measured differential elliptic flow 
carries over from all charged hadrons to identified pions 
and protons. Data for protons that have sufficient statis- 
tical precision to discriminate between different (?7/s)qgp 
values exist only for the mid-centrality range 20 — 50%. 
In very peripheral collisions (> 60% centrality) VISHNU 
has similar problems with the pion V2(px) as we saw in 
Figs. 2 and 3 for all charged hadrons. (We comment 
on this discrepancy in the discussion in Sec. IV.) In the 
mid-centrality range Fig. 5 shows excellent agreement be- 
tween VISHNU with (?7/s)qgp = 0.08 for MC-Glauber and 
(v/ s )qgp — 0-16 for MC-KLN initial conditions and the 
experimental data, in each of the three resp. four cen- 
trality bins shown. The pion elliptic flow data in the left 
column reveal that for both MC-Glauber and MC-KLN 
initial conditions this agreement breaks down if (t]/s)qqp 
is increased by ^ = 0.08 above the preferred value. For 
protons the calculation of V2(px) is numerically costly 
(the elliptic flow signal and the number of protons per 
event are both small), and we have therefore not done 
any calculations for other than the preferred (t?/s)qgp 
values. However, the proton data for these bins are pre- 
cise enough that they would again reject (i]/s)qgp values 
that differed by more than 1/4.7T from the values shown. 



IV. DISCUSSION 



2 Fig. 2 also shows that in both cases the agreement is destroyed 
when increasing (r?/ s )QGP by j- =0.08. 



The comparisons between theory and data discussed 
above prove that we can extract the QGP shear viscosity 
from the centrality dependence of the pt -integrated V2 
for charged hadrons and then use this value to obtain 
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FIG. 5: (Color online) Same as Fig. 2, but for identified pions (left column) and protons (right column). Solid symbols denote 
measurements of «2{2}/(£p art ) 1,/2 from the STAR experiment [19], solid and dashed lines with open symbols show (u2)/e pa rt 
from VISHNU calculations with (?7/s)qgp = 0.08 and 0.16, respectively, using MC-Glauber (top row) and MC-KLN (bottom row) 
initial conditions. Different symbols denote different collision centralities as indicated. 



a very good overall description of the pr- differential V2- 
This works not only for the sum of all charged hadrons, 
but also for individual identified hadronic species, and it 
carries over from V2 to their p^-spcctra which are nicely 
described over the entire range of collision centralities, 
except perhaps the most peripheral collisions. 

We do not recommend to try to extract (t]/s)qqp 
directly from the p^-diffcrential elliptic flow, for the 
following reasons. The main effect of shear viscosity 
that we exploit when extracting it from experiment is 
that it inhibits the hydrodynamic conversion of spatially 
anisotropic pressure gradients within the collision fire- 
ball into momentum anisotropics. As emphasized by 
Ollitrault [29], Heinz [30], and recently by Teaney [31], 
the pr-integrated elliptic flow 3 of the sum of all hadrons 



3 More precisely: the p|,-weighted elliptic flow Ai = P rj 



is the observable that has the most direct relationship 
with the hydrodynamically generated total momentum 
anisotropy. 4 Hence it is the total charged hadron v-i that 
is controlled by rj/s. How the hydrodynamically gener- 
ated total momentum anisotropy is distributed among 
the different hadron species and in pr depends on the 
chemical composition and pt distributions of the hadrons 
[32]. 

The correct theoretical description of the differential 
elliptic flow V2(pt) of individual identified hadron species 
thus depends on the accurate reproduction of their yields 
and pt spectra which show much stronger sensitivities 



/ p —p \ 

rather than V2 = (cos(2</) p )) = ( % | ). 

4 Replacing "all hadrons" by "all charged hadrons" is fine because 
of approximate symmetry between positive, negative and un- 
charged hadrons in ultrarelativistic collisions which generate al- 
most baryon-free fireballs. 
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to details of the hydrodynamic simulation (such as ini- 
tial conditions and shape of the initial density profiles) 
than the total momentum anisotropy itself. For exam- 
ple, in a purely hydrodynamic approach with Cooper- 
Frye freeze-out, even at extremely high collision energies 
where the total momentum anisotropy has time to fully 
saturate before freeze-out, lower freeze-out temperatures 
will lead to more radial flow; this affects the slope of the 
single-particle spectra, causing a concomitant change in 
the slope of V2(pr) which is solely controlled by the fact 
that, after integration over pr, the same total charged 
hadron elliptic flow must be reproduced as for a higher 
freeze-out temperature. In a hybrid approach such as 
ours, the Cooper-Frye procedure used to convert the hy- 
drodynamic output into particle distributions involves a 
so-called "5f correction" [33] describing the deviation 
from local equilibrium on the conversion hypersurface; 
its form is presently not precisely known [34, 35]. For a 
given total charged hadron V2 , different parametrizations 
for Sf lead to different shapes of identified hadron spec- 
tra and V2(pt)- Bulk viscosity has very little effect on 
the total momentum anisotropy (and thus on the total 
charged hadron V2) but affects the radial flow and hence 
the slopes of pr spectra and V2{pt) [35, 36]. These inter- 
dependencies between the hadron pr-spectra and their 
p^-dependent elliptic flow make it hazardous to extract 
(v/ s )qgp from V2{pt)- With such an approach it is 
rather difficult to arrive at a uniformly good description 
of all soft hadron characteristics, and one easily ends up 
with different (t?/s)qgp values extracted from the ellip- 
tic flow of different hadron species or from collisions at 
different centralities. 

When using the pT-integrated charged particle elliptic 
flow to extract (?7/s)qgp one must, however, pay atten- 
tion to the fact that the measured elliptic flow fluctu- 
ates from event to event and may be contaminated by 
non-flow contributions. This was emphasized in Ref. [1] 
where we therefore used elliptic flow data that had been 
corrected for non-flow effects and event-by-event fluctu- 
ations. In the remainder of this article we elaborate on 
how the p^-integrated charged hadron elliptic flow (V2) 
from the dynamical model VISHNU, calculated with the 
( t ?/ s )qgp values extracted in [I], compares directly with 
various experimental measurements that have not been 
corrected for fluctuation and non-flow effects. The trends 
exposed in this comparison provide useful insights. 

Without the ability of doing event- by- event hydrody- 
namic simulations [24, 37, 38], we can at this moment ac- 
count for event-by-event fluctuations of the initial fireball 
density distribution only on average, in one of two ways: 
Either we recenter and rotate each Monte Carlo event, 
in order to align their major and minor axes, before av- 
eraging the density distributions (this produces an aver- 
age density profile s par t in the "participant plane" , char- 
acterized by its average eccentricity £part ~ (£ pa rt) (see 
footnote 1)), or we superimpose the densities without re- 
centering and rotating (producing a smooth average den- 
sity profile srp in the "reaction plane" , with "standard" 



eccentricity e s = |^ + ^| where (in contrast to footnote 
1) the expectation values in numerator and denominator 
are taken with srp). Both of these methods incorporate 
(in different ways) the effect of event-by-event fluctua- 
tions of the shape and orientation of the collision fireball 
on the average initial eccentricity, but do not dynami- 
cally propagate event-by-event fluctuations of the value 
of this eccentricity. As a result, the hydrodynamic evo- 
lution produces a non-fluctuating elliptic flow, and while 
the UrQMD afterburner produces event-by-event V2 fluctu- 
ations, they are only due to finite number statistics and 
not related to event-by-event fluctuations of the initial 
eccentricity e. 

Various experimental techniques measure different va- 
rieties of V2 which are affected in different ways by event- 
by-event flow fluctuations (driven by event-by-event vari- 
ations of the initial eccentricity) and non-flow effects. Hy- 
drodynamic simulations with smooth (non-fluctuating) 
initial conditions indicate a linear relationship V2 oc e for 
not too large eccentricity [23]. 5 If v 2 oc e, the probability 
distribution of the final V2 is directly related t o that of the 
initial e. For example, a measurement of \J would 
yield values that are proportional to yj (e 2 ), with the 
same proportionality constant as between (V2) and (e) 
[39]. Unfortunately, quantities like e pa rt{2} = (Epart) 1 ^ 2 

and £ part {4}= [2(£ 2 art ) 2 - (£ 4 art )] 1/4 that control the 
fluctuation contributions to «2{2} and i>2{4} [20, 21], 
have a different ccntrality dependence than the average 
eccentricities e par t and £ s that characterize our hydrody- 
namic initial conditions. For this reason it has been sug- 
gested in [39] to make comparisons between theory and 
experiment only with appropriately normalized elliptic 
flows. In the absence of non-flow contributions and for a 
linear mapping between e and V2 in each event, one has 
(«2)/(e P art) = «2{2}/e pa rt{2} = t>2{4}/£p ar t{4} where 
the first ratio can be calculated in a single-shot hydrody- 
namic evolution of an average initial profile whereas the 
other two ratios can be measured if the initial eccentricity 
and its fluctuations are known from a model. 

For Gaussian fluctuations, -P(e par t)^ cxp [— ^ £p 2^2 ^ ] > 

1 /2 

it is easy to show [13] that e par t{2} = [(£ pa rt) 2 +c 2 ] re- 
ceives a positive contribution from fluctuations whereas 

£ P art{4} = [(e pa rt) 4 - 2cr 2 (£ part ) 2 - er 4 ] 1/4 is reduced rel- 
ative to (ffpart)- (In fact, we can write in general 

£part{4}= [(e 2 art ) 2 -((4art) ~ (£part) 2 )] 1/4 in terms of 

the difference between two positive definite quantities 



J Event-by-event hydrodynamic simulations with fluctuating ini- 
tial conditions show that flow anisotropy and eccentricity coef- 
ficients of different harmonic order do not completely decouple 
from each other [24, 38], and in both single-shot and event-by- 
event hydrodynamic simulations we have seen evidence [24] that 
in very central collisions the (small) final elliptic flow is affected 
by several harmonic eccentricity coefnents in the initial state 
(which, although all small, are of similar order of magnitude), 
resulting in a nonlinear dependence of V2 on e for small V2- 
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FIG. 6: (Color online) Integrated charged hadron elliptic flow as a function of collision centrality from the PHENIX [44] and 
STAR [19, 43] experiments are compared with VISHNU calculations using participant plane (PP) averaged (a) and reaction 
plane (RP) averaged (b) initial conditions from the MC-KLN and MC-Glauber models and {t]/s)qgp values as indicated. In 
the STAR data and the calculations vi was integrated over the range 0.15GeV/c < pr < 2GeV/c; the PHENIX data were 
integrated over 0.2 GeV/c < pr < 8 GeV/c. 



which, for non-Gaussian distributions, can become neg- 
ative, in which case e par t{4} is not defined.) It has 
been observed in [40] that, for models where e par t shows 
Bessel-Gaussian fluctuations [40, 41], e pa rt{4} agrees ex- 
actly with the reaction plane eccentricity (erp) ~e s , and 
in [21] that i>2{4} is insensitive to two-particle non-flow 
contributions. For these reasons, the authors of [42] 
used hybrid model simulations with reaction-plane aver- 
aged initial conditions for direct comparison with RHIC 
Au+ Au and recent LHC Pb+Pb data [43] . The validity 
of the assumption of (Bcssel)-Gaussian eccentricity and 
flow fluctuations has been challenged in [22] but was re- 
cently validaded for the MC-Glauber and MC-KLN mod- 
els for Au+Au collisions at centralities of up to about 
40% [24]; for more peripheral collisions, the assumption 
breaks down. We here compare results obtained from 
both participant-plane and reaction-plane averaged ini- 
tial conditions with py-integrated i>2{4} data. 

In Fig. 6 we compare VISHNU results with STAR and 
PHENIX data for the integrated charged hadron elliptic 
flow as function of collision centrality. The STAR and 
PHENIX data are integrated over slightly different pr 
ranges; correcting the PHENIX data for the somewhat 
smaller lower pr cutoff used by STAR and in the calcu- 
lations would move them slightly down. In the absence of 
non-flow contributions and the limit of small fluctuations, 
v 2 {EP} « V{v 2 ) 2 +rf and v 2 {4} ps ^/{v 2 ) 2 -a 2 v . Non-flow 
effects would push u 2 {EP} further up but leave v 2 {4} un- 



changed. Calculations with reaction-plane averaged ini- 
tial conditions of eccentricity e part (Fig. 6a) should thus 
fall between i>2{EP} and i>2{4}, perhaps a bit closer to 
«2{4} if the «2{EP} data are affected by non-flow. The 
STAR v 2 {EP] data lie above those of PHENIX, consis- 
tent with the expectation that the PHENIX data should 
have less non-flow contributions (if any at all), but some 
of the difference between the data set (especially at small 
centralities) may also originate from a shift in the central- 
ity definition [26] . Except for the most peripheral central- 
ities, our calculations lie above the PHENIX and roughly 
on the STAR v 2 {EP} data and overpredict the STAR 
v 2 {4} data. This indicates that the chosen (?7/s)qgp val- 
ues (0.08 for MC-Glauber and 0.16 for MC-KLN initial 
conditions) are slightly too small (but not by much, as 
seen by the fact that increasing (?7/s)qgp by 0.08 leads to 
a strong underprcdiction of all data sets), and that the 
slightly larger values of 0.10 for MC-Glauber and 0.20 
for MC-KLN extracted in [1] from fluctuation-corrected 
v 2 data would give better agreement here, too. Fig. 6b 
shows that with reaction-plane averaged initial condi- 
tions the VISHNU results agree very well with the i>2{4} 
data, supporting the argument [39, 40] that (£rp)«£ s 
provides a good substitute for e par t{4}. Again, using the 
slightly larger (?7/s)qgp values from [1] would further im- 
prove the agreement. 

We note the inability of VISHNU to describe the ellip- 
tic flow in the most peripheral collisions where the ex- 
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pcrimentally measured values remain large whereas the 
theoretical predictions decrease rapidly with increasing 
impact parameter. This drop is related to the decreas- 
ing lifetime of the fireball (which is even shorter for runs 
with MC-KLN initial conditions than for Glauber pro- 
files, due to the sharper edges of the MC-KLN profiles 
which lead to faster radial acceleration) . Shorter lifetimes 
leave less time for generating elliptic flow in the fluid dy- 
namic QGP stage, and the highly dissipative hadronic 
stage cannot compensate for this. Calculations with an 
ideal hydro+cascade model that use a different hadronic 
rescattering algorithm (JAM instead of UrQMD) appear to 
share this feature [42]. The fact that the i>2{4} data do 
not show this decrease indicates that non-flow contribu- 
tions (which are not included in the model) are not to 
blame. We do not know how to obtain larger v 2 values 
from the model at large impact parameters. 



V. SUMMARY AND CONCLUSIONS 

We have utilized the recently developed hybrid code 
VISHNU, which couples a relativistic viscous fluid dynam- 
ical description of the quark-gluon plasma (QGP) with 
a microscopic Boltzmann cascade for the late hadronic 
rescattering stage, to calculate charged and identified 
hadron spectra and elliptic flow measured in 200 A GeV 
Au+Au collisions at the Relativistic Heavy-Ion Collider 
(RHIC). We find that, after suitable readjustments of 
initial conditions, the pt spectra of identified hadrons 
(pions and protons) are rather insensitive to the choice 
of the value of the specific shear viscosity rj/s, whereas 
the eccentricity-scaled elliptic flow i^/e shows strong sen- 
sitivity to rj/s. Using initial conditions that incorporate 
event- by- event fluctuations in the initial shape and ori- 
entation of the collision fireball and values (j?/s)qgp for 
the specific shear viscosity of the quark-gluon plasma 
that were recently extracted from the measured central- 
ity dependence of the eccentricity-scaled, p^-intcgratcd 
charged hadron elliptic flow W2,ch/e [1], we were able 
to obtain universally good agreement between theory 
and experiment for the px-spectra and differential ellip- 
tic flow V2(pt) for both pions and protons at all colli- 
sion centralities. Our analysis validates the constraints 
on r)/s reported in our previous work, namely that the 
QGP shear viscosity for T c < T < 2T C lies within the range 
1 < 4tt(t]/s)qqp < 2.5, with the remaining uncertainty 
dominated by insufficient theoretical control over the ini- 
tial source eccentricity s. 
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Appendix A 

In this Appendix we add a few aspects that, due 
to space limitations, were left out from the discus- 
sion in Rcf. [1] of the almost universal dependence of 
the eccentricity-scaled elliptic flow v^js on the charged 
hadron multiplicity density (l/S)dN c \ l /dy on which our 
extraction of the QGP shear viscosity (t]/s)qgp from 
RHIC data was based. Specifically, we show that the 
universality of this dependence (i.e. the feature that it 
only depends on (??/s)qgp but not on any details of the 
initial conditions for the hydrodynamic evolution which 
affect the initial eccentricity e and transverse area S of 
the expanding fireball) holds not only for the participant- 
plane averaged fluctuating initial profiles used in [1] but 
also for the reaction-plane averaged profiles used here in 
Fig. 6b. Furthermore, it is insensitive to the smearing 
area a s used in the MC-Glauber model of Ref. [10] that 
describes the width of the transverse distribution of mat- 
ter created in each nucleon-nucleon collision. However, 
the source eccentriciy e itself depends on this smearing 
area, and hence the (?7/s)qgp value extracted by com- 
paring the universal theoretical i^/e vs. (l/S)dN c h/dy 
curves for different r)/s with a given set of experimental 
V2 vs. dN c h/dy data also depends on this parameter. 

In the traditional MC-Glauber model, one samples the 
positions of nucleons according to the nuclear density dis- 
tributions of the two colliding nuclei and calculates the 
participant eccentricity from the transverse positions of 
the wounded nucleons and/or binary collision points, de- 
scribed by J-functions in the transverse plane. Hirano 
and Nara [10] pointed out that, since the measured nu- 
clear density distribution represents a folding of the dis- 
tribution of nuclcon centers with the finite size of each 
nucleon, the nuclear distribution used for sampling the 
positions of the nucleon centers must be described by dif- 
ferent Woods-Saxon parameters than the measured nu- 
clear density. Following Ref. [10], we therefore use for 
the distribution of nucleon centers in a Au nucleus a 
Woods-Saxon distribution with radius Rau = 6.42 fm and 
surface thickness c?a u = 0.44 fm (instead of the frequently 
used parameters i?£° as - = 6.38 fm and d A ™ as - = 0.535 fm 
that describe the measured nuclear density distribution 
of Au) . We then distribute the entropy of particles emit- 
ted by a wounded nucleon or created in a binary nucleon- 
nucleon collision homogeneously in a cylinder of radius 
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r s = \J a S j-K (where a s is the so-called "transverse smear- 
ing area"), centered at the position of the wounded nu- 
cleon or the collision point and aligned with the beam 
direction. (The same procedure was used in Ref. [22] 
without, however, first correcting the Woods-Saxon pa- 
rameters of the distribution of nucleon centers for the 
finite nucleon size.) 

The value of the smearing area a s is not known a pri- 
ori since it depends on unknown aspects of pre-thermal 
decoherence and entropy production processes. Theoret- 
ically, it is limited from below by the uncertainty prin- 
ciple which does not permit localization of the produc- 
tion points of secondary particles with average trans- 
verse momentum (j>t) to an average distance r s < 1/ (pr) 
from the classical collision point. We explore the choices 
(T s = 42 mb [10, 22] and 4.2 mb. The smaller value is an 
approximation to pointlike secondary particle production 
(our code for calculating the participant-plane averaged 
initial density requires a non- vanishing a s ) ; theoretically, 
it is disfavored by the above uncertainty argument. 

Figure 7 shows the eccentricity (panel (a)) and trans- 
verse area (panel (b)) of the participant-plane averaged 
initial entropy density distribution as a function of the 
number iV part of participant ("wounded") nuclcons in 
Au+Au collisions, for two values of the smearing area 
a s . We see that a smaller smearing area (more pointlikc 
particle production) leads to larger initial fireball eccen- 
tricities e P art and smaller transverse areas S. While the 
effect of varying er s on S is simply an offset, for £ par t 
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FIG. 7: (Color online) Eccentricity £ pa rt (a) and transverse 
area S (b) of the particiant-plane averaged initial entropy 
density distribution from the MC-Glauber model, for smear- 
ing areas a 3 — 42 and 4.2 mb, respectively. (See also Fig. 4 in 
Ref. [22] for comparison.) 



it leads to a change in the slope of its centrality depen- 
dence: changing a s affects £ par t more strongly in periph- 
eral collisions (where the nuclear overlap region is small 
and strongly deformed) than in central ones [22]. 
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FIG. 8: (Color online) (a) Centrality dependence of 
the charged hadron rapidity density per participant pair 
(dN c h/dy) /(N pa , rt / 2). Experimental data are from STAR 
[17] and PHOBOS [45], using dN ch /dy = 1.16 dN ch /drj for 
the PHOBOS data. Theoretical lines are explained in 
the text. (b) Eccentricity-scaled elliptic flow v-zjs as 
function of multiplicity density (l/S)(dN c i 1 /dy), for differ- 
ent values of (??/s)qgp, using participant-plane (PP) and 
reaction-plane (RP) averaged initial entropy density profiles 
from the MC-KLN and MC-Glauber models, normalized to 
(dNch/dy)^* — 810 in the most central (0—5%) Au+Au col- 
lisions. Insensitivity to the smearing area <r 3 in the MC- 
Glauber models is also shown. 

Figure 8 illustrates that changing the smearing area 
in the MC-Glauber model has no effect on the central- 
ity depenence of the produced charged hadron multi- 
plicity (panel (a)) nor on the universality of v^fe vs. 
(l/S)dN c h/dy (panel (b)). In Fig. 8b one sees (see solid 
brown line without symbols for (^/s)qgp = 0.08) that, 
for smaller <7 s = 4.2mb, the reduced fireball area shown 
in Fig. 7b shifts the entire curve towards the right. The 
shift is, however, not horizontal but rather diagonal such 
that, where they overlap, the shifted curve lies on top 
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of the line for the larger value a s = 42 mb. The upward 
component of the shift of the a s = 4.2 mb line arises from 
an increased QGP lifetime, due to the larger initial en- 
tropy density resulting from the smaller initial area 5*; a 
longer QGP lifetime in turn results in a larger momen- 
tum anisotropy at the beginning of the hadronic stage, 
since at RHIC energies the QGP never lives long enough 
for the fireball eccentricity to completely decay before 
hadronization. Larger momentum anisotropy at the be- 
ginning of the hadronic rescattering stage leads to more 
elliptic flow for the finally emitted hadrons. 

In Figure 8 we also show curves obtained from VISHNU 
using reaction-plane averaged (RP) initial profiles instead 
of participant-plane averaged (PP) ones. One sees that 
the centrality-dependence of the final charged multiplic- 
ity per participant (panel (a)) and the dependence of 
V2/S on the multiplicity density (l/S)dN c h/dy are insen- 
sitive to how we average the fluctuating initial profiles 
from the MC-KLN and MC-Glaubcr models when con- 
structing the smooth initial entropy density profile for the 
hydrodynamic evolution. Due to numerical cost we here 
show only curves for shear viscosity (t7/s)qgp = 0.16 for 
both MC-Glauber and MC-KLN models. We have made 
spot checks to convince ourselves that the scaling shown 
in Figure 8 also works for other choices of (t7/s)qgp- 

Figure 9 is a modified version of Fig. 2b in [1] which 
was used to extract the preferred value of (??/s)qgp from 
experimental data by comparing them with VISHNU cal- 
culations using participant-plane averaged initial condi- 
tions from the MC-Glauber model. We remind the reader 
that, even though theoretically the dependence of Vz/e 
on (1/ ' S)dN c h/dy is universal (at least at a fixed colli- 
sion energy [14]) and depends only on a single parame- 
ter, (j]/s)qqp, but not on the initial profile, we do not 
know the correct initial profile that drives the generation 
of elliptic flow in the actual experiments. Since different 
initial conditions have different eccentricities, the same 
set of experimental 1*2 and dN^/dy data yields differ- 
ent Vije and (l/S)dN ( ±/dy when normalized by e and S 
from different initial state models, resulting in different 
extracted values for (?7/s)qgp from a comparison with 
the universal theory curves. The larger eccentricities and 
smaller overlap areas resulting from a MC-Glauber ini- 
tialization with reduced smearing area a s = 4.2 mb lead to 



smaller w 2 /e and larger (l/S)dN c h/dy values (green cir- 
cled data in Fig. 9) than for the MC-Glauber model with 
standard smearing (<7 S — 42 mb, black squares). This re- 
sults in a larger preferred value (t?/s)qgp (closer to 0.16 
than the value of 0.08 we obtained when postulating MC- 
Glaubcr initial conditions with standard smearing) . Fur- 
thermore, normalization of the experimental data with 
MC-Glauber (e, S) values for reduced smearing changes 
the slope of the dependence of v 2 /e on (l/S)dN c h/dy, 
and it no longer agrees with the slope of the universal 
theoretical curves. We conclude that the comparison of 
experimental data with VISHNU results disfavors the hy- 
pothesis that the experimentally measured elliptic flow 
is generated by initial conditions that can be described 
by a MC-Glauber model with almost pointlike secondary 
particle production. This conclusion aligns nicely with 
the theoretical prejudice against such a model on the 
basis that it would violate the uncertainty principle, as 
discussed above. 
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FIG. 9: (Color online) Comparison of the universal theoret- 
ical curves for V2[rj/s]/e vs. (l/S)(dN c h/dy) from the MC- 
Glauber model with a s = 42 mb [1] with experimental data for 
(1)2) [21], normalized by the eccentricity (e par t) and transverse 
fireball area S of the initial profile from the participant-plane 
averaged MC-Glauber model with standard (a 3 — 42 mb) and 
reduced (a s =4.2mb) smearing areas. 
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